subroutine rks3
use var
implicit none

    integer i,m,j
    real(8) du,dv
    real(8) Qmean_n_k(1:gridsx,1:gridsy,5)
    real(8) UU(5,1+Klx:gridsx+krx,1+Kly:gridsy+kry)
    real(8) Ufl(5,0:gridsx,0:gridsy),Ufr(5,0:gridsx,0:gridsy)
    real(8) Ugl(5,0:gridsx,0:gridsy),Ugr(5,0:gridsx,0:gridsy)
        
    call comput_Qx_p_Qx_m
    call comput_Qy_p_Qy_m  
    !UU(1,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,1)
    !UU(2,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,2)
    !UU(3,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,3)
    !UU(4,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,4)
    !UU(5,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,5)
    !call wenojs3(gridsx,gridsy,UU,Ufl,Ufr,Ugl,Ugr)
    !Qx_m(0:gridsx,0:gridsy,1)=Ufl(1,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,2)=Ufl(2,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,3)=Ufl(3,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,4)=Ufl(4,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,5)=Ufl(5,0:gridsx,0:gridsy)
    !
    !Qx_p(0:gridsx,0:gridsy,1)=Ufr(1,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,2)=Ufr(2,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,3)=Ufr(3,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,4)=Ufr(4,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,5)=Ufr(5,0:gridsx,0:gridsy)
    !
    !Qy_m(0:gridsx,0:gridsy,1)=Ugl(1,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,2)=Ugl(2,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,3)=Ugl(3,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,4)=Ugl(4,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,5)=Ugl(5,0:gridsx,0:gridsy)
    !
    !
    !Qy_p(0:gridsx,0:gridsy,1)=Ugr(1,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,2)=Ugr(2,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,3)=Ugr(3,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,4)=Ugr(4,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,5)=Ugr(5,0:gridsx,0:gridsy)
    call Glf(gridsx,gridsy)
    call rhs(Qmean,Qmean_n_k)
    do i=1,gridsx
        do j=1,gridsy
            Qmean(i,j,1)=Qmean_n(i,j,1)-dt*Qmean_n_k(i,j,1)
            
            Qmean(i,j,2)=Qmean_n(i,j,2)-dt*Qmean_n_k(i,j,2)
            
            Qmean(i,j,3)=Qmean_n(i,j,3)-dt*Qmean_n_k(i,j,3)
            
            Qmean(i,j,4)=Qmean_n(i,j,4)-dt*Qmean_n_k(i,j,4)
            !do m=1,4
            !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
            !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
            !    Qmean(i,j,m)=Qmean_n(i,j,m)+dt*du+dt*dv
            !end do            
            Qmean(i,j,5)=Qmean_n(i,j,5)-dt*Qmean_n_k(i,j,5)        
        enddo         
    end do
    call comput_original_var
    call boundary
    
    !k=2
    call comput_Qx_p_Qx_m
    call comput_Qy_p_Qy_m    
    !UU(1,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,1)
    !UU(2,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,2)
    !UU(3,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,3)
    !UU(4,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,4)
    !UU(5,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,5)
    !call wenojs3(gridsx,gridsy,UU,Ufl,Ufr,Ugl,Ugr)
    !Qx_m(0:gridsx,0:gridsy,1)=Ufl(1,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,2)=Ufl(2,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,3)=Ufl(3,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,4)=Ufl(4,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,5)=Ufl(5,0:gridsx,0:gridsy)
    !
    !Qx_p(0:gridsx,0:gridsy,1)=Ufr(1,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,2)=Ufr(2,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,3)=Ufr(3,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,4)=Ufr(4,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,5)=Ufr(5,0:gridsx,0:gridsy)
    !
    !Qy_m(0:gridsx,0:gridsy,1)=Ugl(1,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,2)=Ugl(2,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,3)=Ugl(3,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,4)=Ugl(4,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,5)=Ugl(5,0:gridsx,0:gridsy)
    !
    !
    !Qy_p(0:gridsx,0:gridsy,1)=Ugr(1,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,2)=Ugr(2,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,3)=Ugr(3,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,4)=Ugr(4,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,5)=Ugr(5,0:gridsx,0:gridsy)
    call Glf(gridsx,gridsy)
    call rhs(Qmean,Qmean_n_k)
    do j=1,gridsy
        do i=1,gridsx
            Qmean(i,j,1)=3.d0/4.d0*Qmean_n(i,j,1)+1.d0/4.d0*(Qmean(i,j,1)-dt*Qmean_n_k(i,j,1))
            !write(*,'(I4,I4,2X,100(ES12.4,1X))') i,j,Qmean(i,j,1),Qmean_n_k(i,j,1)
            
            Qmean(i,j,2)=3.d0/4.d0*Qmean_n(i,j,2)+1.d0/4.d0*(Qmean(i,j,2)-dt*Qmean_n_k(i,j,2))
            
            Qmean(i,j,3)=3.d0/4.d0*Qmean_n(i,j,3)+1.d0/4.d0*(Qmean(i,j,3)-dt*Qmean_n_k(i,j,3))
            
            Qmean(i,j,4)=3.d0/4.d0*Qmean_n(i,j,4)+1.d0/4.d0*(Qmean(i,j,4)-dt*Qmean_n_k(i,j,4))
            !do m=1,4
            !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
            !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
            !    Qmean(i,j,m)=3.d0/4.d0*Qmean_n(i,j,m)+1.d0/4.d0*(Qmean(i,j,m)+dt*du+dt*dv)
            !end do
            Qmean(i,j,5)=3.d0/4.d0*Qmean_n(i,j,5)+1.d0/4.d0*(Qmean(i,j,5)-dt*Qmean_n_k(i,j,5))       
        enddo         
    end do
    call comput_original_var
    call boundary
    
    !k=3
    call comput_Qx_p_Qx_m
    call comput_Qy_p_Qy_m   
    !UU(1,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,1)
    !UU(2,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,2)
    !UU(3,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,3)
    !UU(4,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,4)
    !UU(5,1+Klx:gridsx+krx,1+Kly:gridsy+kry)=Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,5)
    !call wenojs3(gridsx,gridsy,UU,Ufl,Ufr,Ugl,Ugr)
    !Qx_m(0:gridsx,0:gridsy,1)=Ufl(1,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,2)=Ufl(2,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,3)=Ufl(3,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,4)=Ufl(4,0:gridsx,0:gridsy)
    !Qx_m(0:gridsx,0:gridsy,5)=Ufl(5,0:gridsx,0:gridsy)
    !
    !Qx_p(0:gridsx,0:gridsy,1)=Ufr(1,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,2)=Ufr(2,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,3)=Ufr(3,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,4)=Ufr(4,0:gridsx,0:gridsy)
    !Qx_p(0:gridsx,0:gridsy,5)=Ufr(5,0:gridsx,0:gridsy)
    !
    !Qy_m(0:gridsx,0:gridsy,1)=Ugl(1,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,2)=Ugl(2,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,3)=Ugl(3,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,4)=Ugl(4,0:gridsx,0:gridsy)
    !Qy_m(0:gridsx,0:gridsy,5)=Ugl(5,0:gridsx,0:gridsy)
    !
    !
    !Qy_p(0:gridsx,0:gridsy,1)=Ugr(1,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,2)=Ugr(2,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,3)=Ugr(3,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,4)=Ugr(4,0:gridsx,0:gridsy)
    !Qy_p(0:gridsx,0:gridsy,5)=Ugr(5,0:gridsx,0:gridsy)
    call Glf(gridsx,gridsy)
    call rhs(Qmean,Qmean_n_k)
    do j=1,gridsy
        do i=1,gridsx
            Qmean(i,j,1)=1.d0/3.d0*Qmean_n(i,j,1)+2.d0/3.d0*(Qmean(i,j,1)-dt*Qmean_n_k(i,j,1))
            
            Qmean(i,j,2)=1.d0/3.d0*Qmean_n(i,j,2)+2.d0/3.d0*(Qmean(i,j,2)-dt*Qmean_n_k(i,j,2))
            
            Qmean(i,j,3)=1.d0/3.d0*Qmean_n(i,j,3)+2.d0/3.d0*(Qmean(i,j,3)-dt*Qmean_n_k(i,j,3))
            
            Qmean(i,j,4)=1.d0/3.d0*Qmean_n(i,j,4)+2.d0/3.d0*(Qmean(i,j,4)-dt*Qmean_n_k(i,j,4))
            !do m=1,4
            !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
            !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
            !    Qmean(i,j,m)=1.d0/3.d0*Qmean_n(i,j,m)+2.d0/3.d0*(Qmean(i,j,m)+dt*du+dt*dv)
            !end do
            Qmean(i,j,5)=1.d0/3.d0*Qmean_n(i,j,5)+2.d0/3.d0*(Qmean(i,j,5)-dt*Qmean_n_k(i,j,5))     
        enddo         
    end do
    call comput_original_var
    call boundary
    
    
    
    

end subroutine 